//Numerical methods are inherited from this class
//
//Joonas Koskinen
//Pro Gradu Työ
//Jyväskylän Yliopisto

#include "cNum.h"

Num::Num() : matrixMethod(false), sparse(false), normFactor(1.0) {
	/*ivaluey=1e-12;
	ivaluety=1e-12;
	ivaluefy=1e-12;
	ivaluetdy=0.0;
	ivaluefdy=0.0;*/

	
	ivaluey=0.00;
	ivaluety=0.00;
	ivaluefy=0.00;
	ivaluetdy=1.0;
	ivaluefdy=1.0;
        
	f_e_tol = 0.1;
}
Num::Num(double ivy, double ivdy) : matrixMethod(false), sparse(false){
	ivaluey=ivy;
	ivaluetdy=ivdy;
	ivaluefdy=ivdy;
}

Num::~Num() {}

//! Return the name of the solver
string Num::getName() {
	return name;
}

//! Set the name of the solver
void Num::setName(string value) {
	name=value;
}

//! Return the vector x
vector<double> Num::get_x() {
	return x;
}

//! Return the vector y
vector <double> Num::get_y() {
	return y;
}

//! Return the vector dy.
vector <double> Num::get_dy() {
	return dy;
}

//! Return the current step size.
double Num::getStep() {
   // std::cout << "returning step = " << step << endl;
	return step;
}

//! Return the x initial value.
double Num::get_ivaluetx() {
	return ivaluetx;
}

//! Return the fx initial value.
double Num::get_ivaluefx() {
	return ivaluefx;
}

//! Return the y initial value
double Num::get_ivaluey() {
	return ivaluey;
}

//! Return the dy initial value.
double Num::get_ivaluetdy() {
	return ivaluetdy;
}

//! Return the fdy initial value.
double Num::get_ivaluefdy() {
	return ivaluefdy;
}

//! Set step size to value.
void Num::setStep(double value) {
    if(value == 0) {
        cout << "Tried to set step to 0" << endl;
        value = 0.01;
    }
    else {
        step=value;
        nstep=-1.0*value;
    }
}

//! Return the tolerance, f_e_tol.
double Num::get_f_e_tol() {
	return f_e_tol;
}

//! Set the tolerance f_e_tol
void Num::set_f_e_tol(double value) {
	f_e_tol = value;
}

//! Set the initial value for x.
void Num::set_ivaluetx(double value) {
	ivaluetx=value;
}

//! Set the initial value for fx.
void Num::set_ivaluefx(double value) {
	ivaluefx=value;
}

void Num::set_ivaluetx2(double value) {
	ivaluetx2=value;
}

void Num::set_ivaluefx2(double value) {
	ivaluefx2=value;
}

//! Set the initial value for y.
void Num::set_ivaluey(double value) {
	ivaluey=value;
}

//! Set the initial value for ty.
void Num::set_ivaluety(double value) {
	ivaluety=value;
}

//! Set the initial value for fy.
void Num::set_ivaluefy(double value) {
	ivaluefy=value;
}

void Num::set_ivaluety2(double value) {
	ivaluety2=value;
}

void Num::set_ivaluefy2(double value) {
	ivaluefy2=value;
}

//! Set the initial value for tdy.
void Num::set_ivaluetdy(double value) {
	ivaluetdy = value;
}

//! Set the initial value for fdy.
void Num::set_ivaluefdy(double value) {
	ivaluefdy = value;
}

void Num::set_ivaluetdy2(double value) {
	ivaluetdy2 = value;
}

void Num::set_ivaluefdy2(double value) {
	ivaluefdy2 = value;
}

//! Set the first value for fx.
void Num::set_fx(double value) {
	fx.clear();
	fx.push_back(value);
}

//! Set the first value for tx.
void Num::set_tx(double value) {
	tx.clear();
	tx.push_back(value);
}

//! Set the first value for fy.
void Num::set_fy(double value) {
	fy.clear();
	fy.push_back(value);
}

//! Set the first value for ty.
void Num::set_ty(double value) {
	ty.clear();
	ty.push_back(value);
}

//! Set the first value for tdy.
void Num::set_tdy(double value) {
	tdy.clear();
	tdy.push_back(value);
}

//! Set the first value for fdy.
void Num::set_fdy(double value) {
	fdy.clear();
	fdy.push_back(value);
}

void Num::tx_push_back(double value) {
	tx.push_back(value);
}

void Num::fx_push_back(double value) {
	fx.push_back(value);
}

void Num::ty_push_back(double value) {
	ty.push_back(value);
}

void Num::fy_push_back(double value) {
	fy.push_back(value);
}

void Num::tdy_push_back(double value) {
	tdy.push_back(value);
}

void Num::fdy_push_back(double value) {
	fdy.push_back(value);
}

void Num::setStepsize(vector<double> stepsize) {
    this->stepsize = stepsize;
}

vector<double> Num::getStepsize() const {
    return stepsize;
}

void Num::setTdy(vector<double> tdy) {
    this->tdy = tdy;
}

vector<double> Num::getTdy() const {
    return tdy;
}

void Num::setFdy(vector<double> fdy) {
    this->fdy = fdy;
}

vector<double> Num::getFdy() const {
    return fdy;
}

void Num::setTy(vector<double> ty) {
    this->ty = ty;
}

vector<double> Num::getTy() const {
    return ty;
}

void Num::setFy(vector<double> fy) {
    this->fy = fy;
}

vector<double> Num::getFy() const {
    return fy;
}

void Num::setDy(vector<double> dy) {
    this->dy = dy;
}

vector<double> Num::getDy() const {
    return dy;
}

void Num::setFx(vector<double> fx) {
    this->fx = fx;
}

vector<double> Num::getFx() const {
    return fx;
}

void Num::setTx(vector<double> tx) {
    this->tx = tx;
}

vector<double> Num::getTx() const {
    return tx;
}

void Num::setGlueIndex(unsigned int glueIndex) {
    this->glueIndex = glueIndex;
}

unsigned int Num::getGlueIndex() const {
    return glueIndex;
}

void Num::setNormFactor(double normFactor) {
    this->normFactor = normFactor;
}

double Num::getNormFactor() const {
    return normFactor;
}

void Num::setLogD(double logD) {
    this->logD = logD;
}

double Num::getLogD() const {
    return logD;
}

//! Normalize the vector's maximum element absolute value to 1.
vector<double> Num::normalize(vector<double> value) {
	double maximum, max, min;
	max=0;
	min=0;
	vector<double> result;
	max = *max_element(value.begin(), value.end());
	min = *min_element(value.begin(), value.end());

	if (fabs(max)>fabs(min)) {
		maximum=fabs(max);
	}
	else
		maximum = fabs(min);

	for (unsigned int i=0; i<value.size(); i++) {
		result.push_back(value[i]/maximum);
	}

	return result;
}

//! Normalize vector y's integral to unity
void Num::norm_y() {
	 y=normalizeInteg(y, x);
}

//! Normalize partial vectors ty and fy.
void Num::normPartials() {
	this->fy=normalize(this->fy);
	this->ty=normalize(this->ty);
}

//! Combine two vectors vec1 and vec2 and return the union.
vector<double> Num::vecCombine(vector<double> vec1, vector<double> vec2, int index){
	reverse(vec2.begin(),vec2.end());
	vec1.insert(vec1.end(),vec2.begin(),vec2.end());
	return vec1;
}

//! Set the initial values for solver.
void Num::initAttribs(cPot * V, cDY * DY) {
	x.clear();
	y.clear();
	tdy.clear();
	fdy.clear();
	ty.clear();
	fy.clear();
	tx.clear();
	fx.clear();
	tx.push_back(ivaluetx);
	fx.push_back(ivaluefx);
	tdy.push_back(ivaluetdy);
	fdy.push_back(ivaluefdy);
	ty.push_back(ivaluety);
	fy.push_back(ivaluefy);
        if(name=="Numerov method") {
            tx.push_back(ivaluetx+getStep());
            fx.push_back(ivaluefx-getStep());
            //tdy.push_back(ivaluetdy2);
            //fdy.push_back(ivaluefdy2);
            tdy.push_back(exp(-1.0*fabs(tx.back())));
            fdy.push_back(exp(-1.0*fabs(fx.back())));
            //ty.push_back(ivaluety2);
            ty.push_back(exp(-1.0*fabs(tx.back())));
            //fy.push_back(ivaluefy2);
            fy.push_back(exp(-1.0*fabs(fx.back())));
            
         //   std::cout << "initattribs: " << ivaluety2 << " " << ivaluefy2 << " " << ivaluefdy2 << " " << ivaluetdy2 << std::endl;
        }
}

void Num::initAttribsPush(double step, vector<double>& ax, vector<double>& ay, vector<double>& ady) {
    ax.push_back(step + ax.back());
    ay.push_back(ivaluefy);
    ady.push_back(ivaluefdy);
}
//! Returns the normed last element of the array value.
double Num::fastnorm(vector<double> value) {
	double maximum, max, min;

	max = *max_element(value.begin(), value.end());
	min = *min_element(value.begin(), value.end());

	if (fabs(max) > fabs(min)) {
		maximum=max;
	}
	else {
		maximum = fabs(min);
	}
	return value.back()/maximum;
}


//! Method to combine all partial vectors to single vectors.
//We clear the partial vectors at the same time to save memory.
//This becomes an issue when the step size is small.
//Multiply the "right side vector" with a factor eta to match these two vectors at a gluepoint.
void Num::combineAttribVectors() {
	double eta=0.0;
	eta = fabs(ty.back() / fy.back());
        
        //If eta == 1.0 there is no need to multiply, only combine.
        if(eta != 1.0){
            for(unsigned int i=0; i<fy.size(); i++) {
                    fy[i]=eta * fy[i];
            }

            for(unsigned int i=0; i<fdy.size(); i++) {
                fdy[i] = eta*fdy[i];
            }
        }
        x.clear();
        y.clear();
        dy.clear();

        x=vecCombine(tx,fx,0);
	y=vecCombine(ty,fy,0);
	dy=vecCombine(tdy,fdy,0);
}

//Method to smooth the vector after method restart at discontinuities
//vector<int> index contains indexes of last elements in calculated blocks
//vector<double> tvec is a vector containing derivatives of vector vec.
void Num::combineAttribVectors(vector<int> index, vector<double> &vec, vector<double> &tvec) {
    double eta = 0.0;
    
    for(unsigned int i=0; i < (index.size() - 1); i++) {
        eta = fabs(vec[ index[i] ] / vec[ index[i] + 1 ]);
        
        for(unsigned int j=index[i] + 1; j <= index[i+1]; j++) {
            vec[j] = eta*vec[j];
            tvec[j] = eta*tvec[j];
        }
    }
}

void Num::clearPartialVectors() {
        tx.clear();
	fx.clear();
	ty.clear();
	fy.clear();
	tdy.clear();
	fdy.clear();
}

void Num::clearAllVectors() {
    y.clear();
    x.clear();
    dy.clear();
    tx.clear();
    fx.clear();
    ty.clear();
    fy.clear();
    tdy.clear();
    fdy.clear();
}
//! Find zeros from the result vector y.
//for each pair {y[i-1],y[i]} i \in [y[1], y.size()-1]
//we check if the signs are different.  If they are
//we increase the "zeros" - counter
int Num::findZeros() {
	int zeros = 0;
	for (unsigned int i=1; i < y.size(); i++) {
		if (y[i-1]>0 && y[i]<0) {
			zeros++;
		}
		else if (y[i-1]<0 && y[i]>0) {
			zeros++;
		}
	}
	return zeros;
}

//! Normalize the argument vector's riemann integral to unity.
vector<double> Num::normalizeInteg(vector<double> value, vector<double> x) {
	vector<double> integrated;
	/*long*/ double sum=0.0;
     
        for(unsigned int i=1; i<value.size()-1; i++) {
            sum = sum + (x[i] - x[i-1])* ((value[i]*value[i]) + (value[i-1]*value[i-1]));
        }

        sum = sqrt(sum * 0.5);

	for(unsigned int j=0; j<value.size(); j++) {
		integrated.push_back(value[j]/sum);
	}
        normFactor = sum;
	return integrated;
}

//! The method of numerical method.  This method is a general wrapper for all numerical methods.
void Num::method(cDY* DY, cPot * V) {
	//Initialize all attributes
        initAttribs(V,DY);
        double nextDiscont = 0.0;
        int indx = 0;
	//Calculate and set the gluepoint
	V->setGlue(V->calculateGluepoint(DY->get_E()));
        
        
        if(V->discontVectorExists()) {
            indx = V->getDiscontVector().size();
            nextDiscont = V->getDiscontVector()[indx];
        }
        else {
            nextDiscont = V->getGlue();
        }
        
    //    cout << "Sizes: tx: " << tx.size() << " ty: " << ty.size() << " tdy: " << tdy.size() << endl;
   //     cout << "Sizes: fx: " << fx.size() << " fy: " << fy.size() << " fdy: " << fdy.size() << endl;
   //    
        
        //Calculate the steps from the right side to the gluepoint.
		double nstep = -1.0*step;
		//while(fx.back() > V->getGlue() /*&& (fx.back() - V->getGlue()) >= step*/) {
//#pragma omp parallel 
//{
  //              #pragma omp single
    //{
                while(true) {
                   // cout << fx.back() << endl;
                  //  while(/*(fx.back() - V->getGlue()) > step &&*/ fx.back() > nextDiscont) {
			try {
			    //If we are close to the gluepoint and still not past it we adjust the stepsize for one step to get exactly to the gluepoint.
                            //cout << "Distance from glue fx.back() -V->getGlue() = " << fx.back() -V->getGlue() << endl;
//                            cout << "Stepsize step = " << step << endl;
			    if((fx.back()-V->getGlue()) < step && (fx.back() - V->getGlue()) > (step/100.0)) {
//                                cout << "Edellinen askel fx: " << fx.back() << endl;
                                //cout << "Askelpituus h(fx) = " << fx.back() -V->getGlue() << endl;
                                double temp_nstep = -1.0*(fx.back() - V->getGlue());
                                methodStep(temp_nstep, fx, fy, fdy, DY, V);
                              //  cout << "Viimeinen askel fx: " << fx.back() << endl;
                                //We have to be careful and check if we use a variable step method.  We may be far away from the gluepoint
                                if((fx.back() - V->getGlue()) < (step/100.0)) {
                                break;
                                }
                            }
                            else if(fx.back() < V->getGlue()) {
                                fx.pop_back();
                                fy.pop_back();
                                fdy.pop_back();
                                //cout << "fx.back() = " << fx.back() << endl;
                                break;
                            }
			    else {
			        methodStep(nstep, fx, fy, fdy, DY, V);
                             //   cout << "fx.back(): " << fx.back() << " fy.back(): " << fy.back() << endl;
                               // cout << "fx.back() - V->getGlue() = " << fx.back() - V->getGlue() << endl;
			    }
			} catch (cNumericalError &ne) {
				cout << "Exception occurred with message: " << endl;
				cout << ne.what() << endl;
				x.clear();
				y.clear();
				tdy.clear();
				fdy.clear();
				ty.clear();
				fy.clear();
				tx.clear();
				fx.clear();
				throw;
			}
                 //   }
                    indx--;
                 if(V->discontVectorExists()) {
                        nextDiscont = V->getDiscontVector()[indx];
                        //init_attribs_push(nstep,fx, fy, fdy);
                       }   
		}
               //cout << "Glue: " << V->getGlue() << endl;
	//Call method_step as long as we are on the left side of the solution (gluepoint is "in the middle").
                indx = 0;
                if(V->discontVectorExists()) {
                        nextDiscont = V->getDiscontVector()[indx];
                }
                else {
                        nextDiscont = V->getGlue();
                }
   // }//end pragma omp single
   // #pragma omp single
   // {
               // while((tx.back() + step) < V->getGlue() /*&& (V->getGlue() - tx.back()) > step*/) {
                while(true) {
                        try {
			    //If we are close to the gluepoint and still not past it we adjust the stepsize for one step to get exactly to the gluepoint.
			    if ((V->getGlue() - tx.back()) < step && (V->getGlue() - tx.back()) > (step/100.0)) {
			        double temp_step = V->getGlue() - tx.back();
//                                cout << "Edellinen askel tx: " << tx.back() << endl;
//                                cout << "Askelpituus h(tx) = " << temp_step << endl;
                                methodStep(temp_step, tx, ty, tdy, DY, V);
                                //cout << "temp_step = " << temp_step << endl;
                                //cout << "Viimeinen askel tx: " << tx.back() << endl;
			    }
                            else if(tx.back() > V->getGlue()) {
                                tx.pop_back();
                                ty.pop_back();
                                tdy.pop_back();
                               // cout << "tx.back() = " << tx.back() << endl;
                                break;
                            }
			    else {
                                methodStep(step, tx, ty, tdy, DY, V);
                                //cout << "tx.back(): " << tx.back() << endl;
                                //cout << "tx.back() < V->getGlue() = " << tx.back() << " < " <<  V->getGlue() << endl;

                                
			    }
			} catch (cNumericalError & ne) {
				cout << "Exception occurred with message: " << endl;
				cout << ne.what() << endl;
                                clearAllVectors();
				throw;
			}
                   // }
                //    cout << V->getGlue()  << " - " << tx.back() << endl;
                    indx++;
                     if(V->discontVectorExists()) {
                        nextDiscont = V->getDiscontVector()[indx];
                     }
		}
  //  } //end pragma omp single
      // cout << "tx.back() = " << tx.back();
     //  cout << " fx.back() = " << fx.back() << endl;
       
    //   cout << "derivative: " << (tdy.back()/ty.back()) - (fdy.back()/fy.back()) << endl;
       //cout << "tdy.size() = " << tdy.size() << " fdy.size() = " << fdy.size() << endl;
    //Yhdistellään vektorit ensin
                //cout << "DONE! " << tx.size() << "  " << fx.size() << endl;
//}//end pragma   
   combineAttribVectors();
  // cout << "derivative after combine: " << (tdy.back()/ty.back()) - (fdy.back()/fy.back()) << endl;
   y = normalizeInteg(y,x);
   
 
     if (fabs(y[ty.size()-1] + y[ty.size()]) < 0.2 && y[ty.size()-1] * y[ty.size()] < 0.0) {
         //   cout << "flippaan" << endl;
            for(unsigned int indx=tx.size(); indx < y.size(); indx++){
		y[indx]= y[indx] * -1.0;
            }
	}

     glueIndex = ty.size() - 1;
     setLogD(tdy.back()/ty.back() - fdy.back()/fy.back());
    //clear partial vectors konmmentoitu testauksen ajaksi
    clearPartialVectors();
    //Normalize partial vectors integrals to unity
	//tdy = normalize_integ(tdy,tx);
	//fdy = normalize_integ(fdy,fx);

}

void Num::methodMatrix(cDY * DY, cPot * V) {
    //The instruction Map<MatrixXd>(arrayd, n, n) = matrix copies the entries of matrix into arrayd. 
    double step = getStep();
    
    std::vector<double> dummy;
    Eigen::VectorXd eigenVec;
    
    methodStep(step, dummy, dummy, dummy, DY ,V);
    cout << "Solving Eigenvalues..." << endl;
    eSolver.compute(H);
    cout << "Done." << endl;
    //y.clear();
    y.resize(H.rows());
    eigenVec = eSolver.eigenvectors().col(0);
    Eigen::Map<Eigen::VectorXd>(&y[0], y.size()) = eigenVec;
    y = normalizeInteg(y,x);
    DY->set_E(eSolver.eigenvalues()[0]);
}

double Num::getEigenValue(int n) {
    
    return eSolver.eigenvalues()[n];
}

void Num::getEigenVector(int n) {
    Eigen::VectorXd eigenVec;
   // cout << "HEP " << n << endl;
    eigenVec = eSolver.eigenvectors().col(n);
    //cout << "eigenVec " << endl;
    y.clear();
    y.resize(eigenVec.rows());
    //cout << "resize " << eigenVec.rows() << endl;
    Eigen::Map<Eigen::VectorXd>(&y[0], y.size()) = eigenVec;
    //cout << "map " << n << endl;
    y = normalizeInteg(y,x);
    //cout << "done " << n << endl;
}

void Num::findNextEigenfunction(cDY* DY, cPot* V, unsigned break_value) {
    //TODO: Siirrettävä muualle setMinE_step();
    double D0,D1;
    double initialEnergy = DY->get_E();
    setMinE_step(1e-9);
    method(DY,V);
    
    D0 = getLogD();
  //  D0 = tdy.back()/ty.back() - fdy.back()/fy.back();
    DY->set_E(DY->get_E() + E_step);
    method(DY,V);
    D1 = getLogD();
    //D1 = (tdy.back()/ty.back()) - (fdy.back()/fy.back());

    while(fabs(D0) < fabs(D1)) {
        DY->set_E(DY->get_E() + getE_step());
        D0 = D1;
        method(DY,V);
        D1 = getLogD();
        //D1 = (tdy.back()/ty.back()) - (fdy.back()/fy.back());
    }

   while(fabs(getE_step()) > getMinE_step() && DY->get_E() > initialEnergy) {

        while(fabs(D0) > fabs(D1)) {
            DY->set_E(DY->get_E() + getE_step());
            D0 = D1;
            method(DY,V);
            D1 = getLogD();
            //D1 = (tdy.back()/ty.back()) - (fdy.back()/fy.back());
            //cout << "D0 = " << D0 << " D1 = " << D1 << " : fabs(D0) > fabs(D1) = " << (fabs(D0) > fabs(D1)) << endl;
        }

        setE_step(-0.5*(getE_step()));
        DY->set_E(DY->get_E() + getE_step());
        //cout << "E = " << DY->get_E() << endl;
        D0 = D1;
        method(DY,V);
        D1 = getLogD();
        //D1 = (tdy.back()/ty.back()) - (fdy.back()/fy.back());
    }
   // cout.precision(15);
   // cout << "Energy = " << DY->get_E() << " D0 = " << D0 << " D1 = " << D1 << " tx.back() = " <<  tx.back() << " fx.back() = " << fx.back() << endl;
}

void Num::shooting(cDY* DY, cPot* V) {
    initAttribs(V, DY);
 //   double nstep = -1.0 * step;

    while(E_step > 1e-15) {
        E_step = 0.5 * fabs(E_step);
        cout << "E = " << DY->get_E() << "  E_step = " << E_step << endl;
        while(y.back() <= y[0]){
            DY->set_E(DY->get_E() + E_step);
            while(x.back() < V->getUniverse()) {
                methodStep(step, x, y, dy, DY, V);
            }
            cout << "E = " << DY->get_E() << endl;
        }
        E_step = -0.5 * fabs(E_step);
        while(y.back() > y[0]) {
           DY->set_E(DY->get_E() + E_step);
           cout << "E = " << DY->get_E() << endl;
           while(x.back() < V->getUniverse()) {
                methodStep(step, x, y, dy, DY, V);
            }
       }

        if(y.back() == y[0]) {break;}
    }
        
    cout << "y[" << x[0] << "] = " << y[0] << endl;
    cout << "y[" << y.size()-1 << "] = " << y.back() << endl;
    cout << "E = " << DY->get_E() << endl;
    cout << "pot right = " << V->getPotRight() << endl;
}

void Num::setMatrixMethod(bool matrixMethod) {
    this->matrixMethod = matrixMethod;
}

bool Num::isMatrixMethod() const {
    return matrixMethod;
}

void Num::setH(Eigen::MatrixXd H) {
    this->H = H;
}

Eigen::MatrixXd Num::getH() const {
    return H;
}

void Num::setY(vector<double> y) {
    this->y = y;
}

vector<double> Num::getY() const {
    return y;
}

void Num::setX(vector<double> x) {
    this->x = x;
}

vector<double> Num::getX() const {
    return x;
}

void Num::setEigenIndex(int eigenIndex) {
    this->eigenIndex = eigenIndex;
}

int Num::getEigenIndex() const {
    return eigenIndex;
}

void Num::setSparseH(Eigen::SparseMatrix<double> sparseH) {
    this->sparseH = sparseH;
}

Eigen::SparseMatrix<double> Num::getSparseH() const {
    return sparseH;
}

void Num::setSparse(bool sparse) {
    this->sparse = sparse;
}

bool Num::isSparse() const {
    return sparse;
}

//! Check if value1 and value2 have the same sign.
bool Num::isSignSame(double value1, double value2) {
    if(signbit(value1) == 0 && signbit(value2) == 0) {
        return true;
    }
    else if(signbit(value1) != 0 && signbit(value2) != 0) {
        return true;
    }
    else if( ( signbit(value1) != 0 && signbit(value2) == 0 ) || ( signbit(value1) == 0 && signbit(value2) != 0 ) ) {
        return false;
    }
    cout << "Ei saa olla täällä" << endl;
    return false;
}
void Num::setNextNode(unsigned int nextNode) {
    this->nextNode = nextNode;
}

unsigned int Num::getNextNode() const {
    return nextNode;
}
void Num::setMinE_step(double minE_step) {
    this->minE_step = minE_step;
}
double Num::getMinE_step() const {
    return minE_step;
}
void Num::setE_step(double E_step) {
    this->E_step = E_step;
}
double Num::getE_step() const {
    return E_step;
}

double Num::getMaxerror() const
{
    return maxerror;
}

double Num::getMaxstep() const
{
    return maxstep;
}

double Num::getMinerror() const
{
    return minerror;
}

double Num::getMinstep() const
{
    return minstep;
}

void Num::setMaxerror(double maxerror)
{
    this->maxerror = maxerror;
}

void Num::setMaxstep(double maxstep)
{
    this->maxstep = maxstep;
}

void Num::setMinerror(double minerror)
{
    this->minerror = minerror;
}

void Num::setMinstep(double minstep)
{
    this->minstep = minstep;
}

//! Normalize two vectors with the same factor and insert normalized last elements to na and nb.
void Num::normalize_vectors(vector<double> a, vector<double> b, double & na, double & nb) {
	vector<double> max;

	max.push_back(fabs(*max_element(a.begin(), a.end())));
	max.push_back(fabs(*min_element(a.begin(), a.end())));
	max.push_back(fabs(*max_element(b.begin(), b.end())));
	max.push_back(fabs(*min_element(b.begin(), b.end())));

	na = a.back() / *max_element(max.begin(), max.end());
	nb = b.back() / *max_element(max.begin(), max.end());
}

//! Return the i:th element of the vector x.
double Num::get_xi(unsigned int i) {
	return x[i];
}
//! Return the i:th element of the vector y.
double Num::get_yi(unsigned int i) {
	return y[i];
}

//! Return the size of the vector x.
unsigned int Num::sizeof_x(){
	return x.size();
}

//! Return the size of the vector y.
unsigned int Num::sizeof_y(){
	return y.size();
}

bool Num::isResultLimitSet() {
	return resultLimited;
}

void Num::useResultLimit(bool value) {
	resultLimited = value;
}

double Num::getResultLimit() {
	return resultLimit;
}

void Num::setResultLimit(double value) {
        if(y.size() > 0) {
            resultLimit = value;
//            useResultLimit(true);
            limitResultVector();
        }
     //   else {cout << "Resultlimit not set. No result vector." << endl;}
}

unsigned int Num::getLimitIndexBegin() {
	return limitIndexBegin;
}

unsigned int Num::getLimitIndexEnd() {
	return limitIndexEnd;
}

void Num::limitResultVector(void) {
	setLimitIndexBegin(0);
	setLimitIndexEnd(y.size()-1);

         if(isResultLimitSet() == true){
                cout << "reslim is true" << endl;
             for(unsigned int i=0; i < y.size(); i++) {
			if(fabs(y[i]) > getResultLimit()) {
                            setLimitIndexBegin(i);
                            break;
			}

		}
		for(unsigned int i = y.size()-1; i > 0; i--){
			if(fabs(y[i]) > getResultLimit()) {
                            setLimitIndexEnd(i);
                            break;
			}
        }

	}
}

void Num::setLimitIndexBegin(unsigned int value) {
	limitIndexBegin = value;
}

void Num::setLimitIndexEnd(unsigned int value) {
	limitIndexEnd = value;
}
